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We use a self-consistent strong-coupling expansion for the self-energy (perturbation theory in 
the hopping) to describe the nonequilibrium dynamics of strongly correlated lattice fermions. We 
study the three-dimensional homogeneous Fermi-Hubbard model driven by an external electric field 
showing that the damping of the ensuing Bloch oscillations depends on the direction of the field, 
and that for a broad range of field strengths, a long-lived transient prethermalized state emerges. 
This long-lived transient regime implies that thermal equilibrium may be out of reach of the time 
scales accessible in present cold atom experiments, but shows that an interesting new quasi-universal 
transient state exists in nonequilibrium governed by a thermalized kinetic energy but not a thermal- 
ized potential energy. In addition, when the field strength is equal in magnitude to the interaction 
between atoms, the system undergoes a rapid thermalization, characterized by a different quasi- 
universal behavior of the current and spectral function for different values of the hopping. 

PACS numbers: 03.75.Ss, 03.65.Yz, 71.10.Fd 



Introduction. Nearly all experiments with cold atoms 
involve nonequilibrium processes in either the experimen- 
tal probes used for measurements or in the experimental 
setup itself. Of particular interest is the nonequilibrium 
behavior when the interactions are strong enough to cre- 
ate a Mott insulator. In this regime, perturbative meth- 
ods based on an expansion in the interaction strength are 
bound to fail. In this work, we consider a different, com- 
plementary technique, based on a strong-coupling expan- 
sion [3-IH , and demonstrate its application to nonequilib- 
rium systems. The strong-coupling expansion has proved 
to be an accurate and well controlled approximation for 
equilibrium cold atom problems [H, Q (especially for tem- 
peratures larger than the hopping), yielding excellent ap- 
proximate results at a fraction of the computational cost 
of more exact methods. 

In nonequilibrium, it is crucial to account for the 
damping and relaxation effects in a perturbed system. 
Typical random phase approximation (RPA) type cal- 
culations include the lowest-order correction to the self- 
energy and cannot capture these damping effects. So we 
must go to the next order and we include all corrections 
up to the second order in hopping, which allows us to cap- 
ture the damping effects, although at an increased com- 
putational cost. Other related techniques for the strong- 
coupling regime include the hybridization expansion [5| 
used in dynamical mean- field theory (DMFT), but our 
approach will work in any dimension. 

To illustrate this method, we study Bloch oscillations 
in a Mott insulating cold atom system in an artificial elec- 
tric field. Similar work has been done within the context 
of DMFT Here we show that a prethermalized 



state persists for a long time after the transient exci- 
tation of the system, and that it shows quasi-universal 
behavior for different strengths of the field or hopping 
(similar prethermalized states have been studied in other 
contexts [Iol - fl2l ] ) . We also find a very rapid thermaliza- 
tion when the field equals the interaction strength with 
scaling behavior for the current and for the density of 
states for a wide range of hoppings. 

Model. We describe the nonequilibrium cold atom 
system with a single-band Fermi-Hubbard model for a 
homogeneous system with no trap potential: 

H(t) =j2 H i 0) ( t )+ Hihop) ( t ) 

i 

H{t) =Y^U(t)n W - J(t) ]T (e a ^clc ]a + h.c.) , 
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where J(t) is the hopping amplitude between the nearest 
neighbor lattice sites, U (t) is the energy penalty for two 
atoms to occupy the same optical lattice site, A(t) is the 
vector potential fully describing the external "electric" 
field E(t): A(t) = - f* E(t')dt' (which corresponds to 
"pulling" the optical lattice through the atomic cloud in 
a cold atom experiment). c\ a is the creation operator for 
an atom on site i in the hypcrfine ("spin") state a with 
two available spin states: a =t, \.. The occupancy at site 
i in state a is = c\ a ci a and = fi — fj, where fi 
is the position vector of lattice site i. In this model, the 
hopping J(t), interaction U(t), and field E{t) can be time 
dependent. The vector potential is used to describe the 
electric field because it preserves translation invariance. 
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Due to gauge invariance, the same field can be described 
by an appropriate scalar potential, which corresponds to 
a "lattice tilting" by adding a linear potential in a cold 
atom experiment. 

Methodology. We solve the model using the strong- 
coupling expansion [H, 0, 01, which treats the hopping 
amplitude small parameter. The single-particle 

propagation is fully described by the contour-ordered 
Green's function G(t,t') = -i(T c c(t)J (*')) where both 
time arguments lie on the Kadanoff-Baym-Kcldysh con- 
tour (Fig. [TJ. In the atomic limit, the Green's function 
on site i is given by Gf\t,t') = — i(T c a(t)cl(t')) H (o) . 
The first-order correction to the full Green's function is 
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FIG. 1: The Kadanoff-Baym-Keldysh contour describing the 
nonequilibrium evolution of the system in real time t — 
. . . tmax following an initial (equilibrium) thermalization at 
temperature T. Numerical solution requires discretization re- 
sulting in N t = 2N[ e + Nl m finite time steps. 



due to a single hopping to/from a neighboring site: 
GfMt,t')= [ dt 1 Gf\tM)JiAti)Gf\h,t'), (2) 



where the time integral is taken over the full contour. All 
local second-order terms can be expressed as [13| : 



where 



Gi 2) (M') = f dti f dt 2 ^J im (t 2 )J m i(ti) 
Jc Jc m 

0i o) (*o,ti;t 3 ,t3) = ai o) (tb,*i,*2,*3) 
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is a second-order cumulant Green's function and 
Qi 0) (to,ti,t2,t 3 ) = -i(T c c l (io)c 4 (ti)c l t (i 2 )c 4 t (i 3 )) ff (o) is 
the two-particle atomic Green's function. 

Due to the increased computational complexity, the 
terms beyond second order are truncated. However, a 
partial but infinite resummation of terms is possible. For 
a homogeneous system, this leads to a matrix equation 
for the momentum-dependent Green's function: 



G 



k = {[GW]- 1 -5 c J k -V[G]] 1 , 



where Jk = — 2 J(t) J2i=i cos(ki) is the Fourier trans- 
form of the hopping on a c?-dimensional cubic lattice, 
and we have abandoned the time indices in favor of a 
matrix notation. This requires discretization of the con- 
tour, so that the matrix size equals the number of time 
slices on the contour (Nt). We adopt a convention that 
G(t,t) = G > (t,t) for the equal time Green's function, 
which results in following definition of the delta function 
on the discretized contour S c .ij = 6i+ij — 5i,i5j.N t f° r 
i, j = l...N t . Lastly, S[G] is the "second-order cumulant 
self-energy" : 
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which implicitly depends on the local single-particle 
Green's function on neighboring sites used to calculate 

(2) 

G u . The only non-local second-order term (which cor- 
responds to hopping two sites away) does not appear in 
the expression in Eq. ([S]), since it is recovered by resum- 
mation as a product of first-order terms. Finally, we 
require that the propagation on the adjacent site in the 
self-energy functional (Eq. [5]) is described by the same 
full Green's function as on the original site, or, equiva- 
lently, that the local Green's function obtained through 
the momentum summation is the same as the one used 
to calculate the self-energy S[G]. This approximation is 
justified in a homogeneous system and results in a so- 
lution that captures the relaxation effects to a steady 
state. It also imposes a self-consistency condition which 
requires an iterative solution. Remarkably, this method 
also recovers the exact result for [7 = 0, since in that case 
both the second-order cumulant Green's function and the 
self-energy vanish. A diagrammatic interpretation of the 
approximation is shown in Fig. [2] 
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FIG. 2: a) Graphic notation for the full Green's function, 
atomic Green's function, hopping and the second-order local 
Green's function (Eq. [3}. b) Diagrammatic interpretation of 
the self-consistency condition. The resummation in momen- 
tum space requires that the initial propagation in the first and 
second-order terms is governed by the full Green's function. 
In addition, we require that the propagation on the neighbor- 
ing site in the second-order term is also the same as the full 
Green's function. 



A number of measurable quantities can be obtained di- 
rectly from the single-particle Green's function. The mo- 
mentum distribution is calculated as: rik(t) = iGk(t u , t{), 
where t u (t{) denote upper (lower) branches of real time 
part of the Keldysh contour (Fig. [IJ . The kinetic energy 
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and current can be calculated as sums over fc-space 141 ) : 

N k d 



e km {t) 



2J(t) 



^2 X] cos ( fcm - A m (i))n fc (t) (7) 



2J(t) Nk 

3m(t) = -j^ 1 sin (k m - A rn (t))n k (t) . (8) 



The potential energy (and double occupancy) can be ob- 
tained from the equal time derivative of the Green's func- 
tion and the previously calculated kinetic energy: 



dG{t,t') 



at 



nt(*K(*)-|(»t(*) + ^(*)) + 5 



e km (t) 



(9) 



The local retarded Green's function is: 

G R {t,t') = -iQ{t-t'){{c{t),J (*')}) 
= G(ti,t u ) - G{t u ,t[) , 



(10) 



and the spectral function at an average time i Q is 
A ta (w) = -ilm / °° G R (t a + t/2, t a - t/2)e^dt. 

Results. We take the hopping, interaction and elec- 
tric field to be time independent. We investigate a three 
dimensional cubic system of N = 32 3 sites at half filling 
(the total number of atoms is equal to the number of op- 
tical lattice sites). The hopping J/U is taken small so 
that the system is originally in a Mott-insulating state. 

We start with a system initially in thermal equilibrium 
at a temperature T/U = 0.25 (corresponding to 6.5% of 
sites doubly occupied) and turn on a constant field at 
time t = 1. The resulting response of the energy and 
the current is shown in Fig. [3l Sudden application of 
a strong field (E > U) gives rise to Bloch oscillations 
with a period of 2-k/E. The damping or decay of the 
Bloch oscillations depends on the number of axial direc- 
tions for which the field is nonzero. For a field aligned 
with the axial direction of the lattice, the Bloch oscilla- 
tions decay rapidly with a time-scale determined by the 
inverse hopping (J -1 ), whereas for a field aligned with 
the diagonal direction, there is essentially no damping 
within the time scale of the simulation. A strong field 
suppresses the tunnelling of particles along the direction 
of the field, resulting in an effective two-dimensional sys- 
tem for the axial field, and an effective zero-dimensional 
system for a diagonal field [HJ. The amplitude of the 
current is proportional to the hopping parameter J. We 
show results for one initial temperature here. As the 
initial temperature increases, the amplitude of the cur- 
rent decreases. The additional amplitude modulation of 
the current has a dominant frequency of U, that is inde- 
pendent of field, temperature or hopping. For both field 
directions, a strong field results in a small increase of the 
total energy that is almost entirely due to the change in 



the kinetic energy which is much smaller than the energy 
in the long-time thcrmalized state, which will have the in- 
finite temperature limit of the energy. Clearly, the time- 
scale for full thcrmalization of the system far exceeds 
the time of simulation, as expected, since the creation of 
double occupancies involves multiple particle processes 
when the hopping is much smaller than U (due to energy 
conservation). 
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FIG. 3: (color online) Total energy (top) and current (bot- 
tom) for different values of electric field at U/T = 4 and 
U/J = 24. The field is suddenly switched on at t = 1. A 
strong field results in Bloch oscillations in the current and 
the total energy, yet the damping of the oscillations depends 
strongly on the direction of the field. The dotted line shows 
the infinite temperature equilibrium limit for the total energy 
(Etot = U/4), and the dashed lines correspond to equilibrium 
with the initial temperature and no field. The current for 
cases when the field is aligned with the axial direction of the 
lattice is scaled by \/3 for a better comparison with the results 
for a diagonal field. 



Setting the field equal to the interaction (E — U) leads 
to markedly different behavior illustrated by a rapid rise 
of total energy that approaches (or even overshoots) the 
maximum (T — > oo) value for both directions of the field. 
The current shows a peak that is followed by a slow decay 
with a time-scale on the order of J™ 1 . No Bloch oscilla- 
tions are present, but a slight period 2ir/U modulation is 
apparent during the initial rise of the current. Scaling of 
the current and time for different values of hopping re- 
veals the quasi-universal nature of the decay (see Fig. [4}. 
The overall shape of the response depends on the direc- 
tion of the field and the dimensionality of the system. 
Smaller fields (E < U) again lead to Bloch oscillations, 
which become less pronounced with a further decrease of 
the field 0. 

We show the spectral function for the prethermalized 
state in Fig. [5] for a diagonal field of variable strength. 
Even for a small field, the spectral function rapidly col- 
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FIG. 4: (color online) Scaled current for different values 
of hopping for E = [7(1,0,0) and T = J, showing quasi- 
universal behavior. The data change very little with a further 
decrease of the initial temperature. The exponential decay of 
the current has a time-scale of J -1 . The inset shows scaled 
spectra (upper Hubbard band) at an average time t a = 1/J, 
which also show a quasi-universal behavior. 



lapses into sharp peaks near lu/U — ±0.5 with smaller 
sidebands split off the main peaks roughly by plus or 
minus the size of the electric field. For E/U = 0.5, a res- 
onant peak emerges at zero frequency, similar to results 
obtained by a numerical rcnormalization group method 
in fT3] - In the E w U case, the Hubbard bands widen due 
to rapid thermalization. Since the retarded Green's func- 
tion data are restricted to short time intervals, we use a 
method similar to the one introduced by Sandvik Il6| to 
transform data to the frequency domain (see also [13j|). 
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FIG. 5: (color online) Spectral function at an average time 
t a = 32.0 for U/J = 24 and U/T = 4 for a diagonal field. Due 
to the even symmetry of the spectral function for the half filled 
system, only positive energies are shown. The application of 
a field results in a sharp peak near u/U = 0.5 in all cases, 
except when E/U — 1.0. For E/U = 0.5, a zero frequency 
peak emerges. See Ref. [ll|] for details on numerics. 



When the field is aligned with the lattice, the Bloch 
oscillations arc rapidly damped and the Hubbard bands 
remain wide and featureless [13j . except for E/U = 0.5, 
which shows a spectral weight transfer from co/U = ±0.5 
to a zero frequency resonance (l7j . 

The prethermalized state persists for long times (longer 
than the inverse hopping J -1 ), and unfortunately we 
cannot directly measure the associated time-scale which 
is probably exponential in U/J @ or longer, and has 
been sugg ested as related to the slow creation of dou- 
blons In a strong field, all our observations are 

consistent with a dimensional reduction mechanism sug- 
gested in (l5j |. 

Conclusions. Based on a strong-coupling nonequilib- 
rium method, we have established the generic presence 
of a long prethermalized regime in a Mott insulator for 
a wide range of uniform field strengths. While damping 
of Bloch oscillations strongly depends on the direction 
of the field relative to the axial directions of the lattice, 
fast thermalization occurs only when field strength ap- 
proaches the repulsive interaction between atoms. In 
this case, the time scale of the thermalization is de- 
fined by the hopping, and several measurements, includ- 
ing the single-particle spectrum, point to quasi-universal 
behavior, which is independent of hopping or tempera- 
ture (within a certain range). For field strengths different 
from the interaction, the heating is strongly suppressed 
and the system enters a prethermalized regime charac- 
terized by a long lifetime which unfortunately cannot be 
computationally resolved. The striking difference in the 
time scales of thermalization for different field strengths 
are experimentally relevant. The strong-coupling method 
introduced in this work is quite general and can be used 
to study a wide range of noncquilibrium phenomena. 
Despite an inability to reach the lowest temperatures, 
it is especially well suited to study cold atom systems, 
where current experiments have similar limitations. The 
method is also easy to extend for studies of bosonic sys- 
tems. In future work, we plan to generalize it to include 
the trap potential. 
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